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Abstract 



We study the iron dimer by using an accurate ansatz for quantum chemical calculations based on a simple variational wave function, 
defined by a single geminal expanded in molecular orbitals and combined with a real space correlation factor. By means of this 
approach we predict that, contrary to previous expectations, the neutral ground state is 7 A while the ground state of the anion is 8 X~, 
hence explaining in a simple way a long standing controversy in the interpretation of the experiments. Moreover, we characterize 
consistently the states seen in the photoemission spectroscopy by Leopold et al..\M It is shown that the non-dynamical correlations 
included in the geminal expansion are relevant to correctly reproduce the energy ordering of the low-lying spin states. 
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1. Introduction 

The iron dimer is a puzzling molecule. Indeed, explaining 
all the related experimental findings in a consistent theoretical 
frame is very hard. Moreover, this dimer represents a testcase to 
check the ability of a method to capture the physics of the tran- 
sition metal compounds, since it includes all their difficulties, 
namely a strong electron correlation in the nearly half-filled d 
orbitals, and a non trivial ordering of the low-lying energy states 
' which differ by their spin. 

In 1986 Leopold et al. carried out an experiment HI of neg- 
ative ion photoelectron spectroscopy (PES), with the aim of 
studying the low-lying electronic states of Fe2. A sample of 
Ftj is prepared and excited by an incoming photon. The spec- 
trum of Fe2 appears remarkably simple, with only two peaks, 
corresponding to the excitations from the Fe^ ground state to 
those Fe2 states allowed by the selection rules, which implies 
that the total angular momentum of the final state cannot change 
by more than one. Both states reveal the same vibrational fre- 
quency and bond length. Few years later, Leopold[2] argued 
that the simplest explanation of these data is to admit that the 
9 Y, g is the ground state of Fe2, and the 8 E~ is the ground state of 
Fe 2 . This interpretation is based on the hypothesis that the two- 
band system observed in the Fe^ spectrum is due to the detach- 
ment from a 4s-like molecular orbital (MO). Therefore, if one 
supposes that the ground state of Fe~ is 8 E~, its configuration 
turns out to be cftAs) cr* u ' 2 (4s)3d 13 , and the 4s electron detach- 
ment would produce two possible states, with the same orbital 
configuration cr 2 (4s) o~ u ' (4s)3d 13 but with the <r* u {4s) of high 
( 9 l, g ) or low ( 7 2p spin coupled to the remaining 3ii 13 electrons. 
These two states would correspond to the first and second peak 
of the spectrum, respectively, and display the same structural 
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properties, the only difference being the spin coupling. How- 
ever, this interpretation disagrees with electron spin resonance 
experiments, 13|] which fail to observe the Fe2 and therefore sup- 
ports the idea of an orbitally degenerate ground state, like the 
7 A„ configuration, unless the iron dimer features a large zero- 
field split (larger than 8 cm -1 ), thus producing an energy split 
not detectable by the experimental setup. However, such a large 
magnetic split is unusual for that kind of molecules. 

From the theoretical side, some numerical studies based 
mainly on density functional theory (DFT) methods with vari- 
ous functionalsfl4j, |5l [g, |7[| and multi reference configuration in- 
teraction (MRCI) calculations [8] yielded the 7 A U as the ground 
state of Fe2, while more recent MRCI calculations[9] and DFT 
studies with coupled cluster lfioll and +t/ lUlll corrections sup- 
ported the idea that its ground state is 9 Z~. Those methods gave 
a 8 Ag ( 8 X~) ground state for the anion whenever a 7 A„ ( E~) 
ground state was found for the neutral dimer. 

In this Letter, we tackle the study of the iron dimer by means 
of quantum Monte Carlo (QMC) simulations based on the res- 
onating valence bond (RVB) wave function. 1 12l 1 1 311 which is a 
correlated single-determinant ansatz successfully used in previ- 
ous calculations. lfl4l \l5 , HI 17] Here, we use an extension of 
the RVB picture that is based on a MO expansion of the sin- 
glet state in the determinant. By setting the number of MO's to 
a value such that a Jastrow correlated single determinant (SD) 
wave function is recovered for each fragment in the atomization 
limit, we obtain a description of the bond which is remarkably 
accurate. This ansatz has been tested on a set of first-row atom 
dimers, such as F2, N2, Be2, and C2, where our calculations 
yielded results for the binding energy and the equilibrium dis- 
tance very close to the experimental values . 11 1811 much better 
than those obtained in previous calculations based on standard 
Jastrow-SD (JSD) wave functions.|fl9| 
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2. Method 

We consider a RVB wave function for N electrons, where for 
simplicity we take Nf > iVj,, and iVj (N{) is the number of spin 
up (down) electrons. The wavefunction is given by the product 
of a determinantal part and a Jastrow correlation factor. The de- 
terminantal part is the antisymmetrized product of singlet pairs. 
Each pair is described by a symmetric (p(r, f) = <p(? , r) orbital 
function. In order to define a pure spin state with total spin 
S — \N-[ - N[\/2 and maximum spin projection S'"' = S , we 
take Ni singlet pairs and 2S unpaired orbitals <f)j(f) for spin up 
electrons, and apply the antisymmetrization operator J\ to ob- 
tain a consistent fermionic wave function: 



_ i i 



(i) 



where <$>n is the A^-electron wave function and R indicates 
the corresponding 3A^-dimensional vector of coordinates, R = 



f[, ■ ■ ■ , r" N , r[, ■ • • , r" Ni J. The wave function in Eq. ([TJ is also 



J, 

called antisymmetrized geminal product (AGP) and can be 
computed by means of a single determinant (see Ref. 11411 and 
references therein). 

As said, the total spin of <C>^(/?) is definite. We also impose 
all possible symmetries to be satisfied, including angular mo- 
mentum and spatial reflections. 

The Jastrow correlation factor is the other important ingredi- 
ent of the wavefunction. Its generally adopted form reads: 



J(x) - exp 



V «J 



(2) 



where f(f,f) is a function of two electron coordinates.] 15] 
The Jastrow term accounts for the electron-electron repulsion 
and suppresses configurations with overlapping valence bonds, 
which would lead to a too large electron density around an 
atom, with an increase in the total energy. 

As any function of two coordinates, the pairing function (p 
as well as the correlation function / in the Jastrow term can be 
expanded in terms of single particle orbitals. In particular, the 
pairing function reads: 



n-ZS 



(3) 



where n is large enough, and {<pj} is an orthogonal single par- 
ticle basis se^, which reaches its complete basis set limit for 
n — > oo. Notice also that in these notations we assume that the 
2S unpaired orbitals ; correspond to the indexes: n - 2S + 1 < 
; ^ n. 

The single particle orbitals <pj can be conveniently chosen 
as the MO's obtained with a conventional restricted Hartree- 
Fock (RHF) calculation. Indeed the MO basis allows us to write 



'The orbitals <f>i can be MO's expanded in terms of atomic orbitals tp a j 
where a indicates the atomic center and j the type: tf>i(f) = Yia,jX' a jfajif)- 
The coefficients x' a j, as well as the weights Aj, can be used as variational pa- 
rameters defining the geminal in Eq. {3} 



Eq. (0 in a diagonal form equivalent to a more involved matrix 
form when the MO's are developed in an atomic basis set. lfl5ll 
By truncating the expansion in Eq. (f3]l to a number of MO's 
n equal to the number of electron pairs and unpaired orbitals, 
namely n = N^, one recovers the usual RHF theory, because 
the antisymmetrization operator J\ clearly singles out only one 
Slater determinant. Moreover, the MO weights Aj affect only an 
overall prefactor of this Slater determinant, so that their actual 
values are irrelevant in this case. However, the pairing function 
is generally not limited to have only N± non vanishing eigen- 
values Aj. Therefore, the RVB wave function represents a clear 
extension of the RHF theory, not only for the presence of the 
Jastrow factor, which considerably improves the dynamical cor- 
relations, but mainly because its determinantal part goes beyond 
RHF when n > Nf, by including also non-dynamical correla- 
tions. Quite generally, a gain in energy and a more accurate 
calculation are expected whenever n > Nf. 

In this Letter, we use in all calculations a number n* of MO's 
that is enough to have a fully symmetric state that connects the 
compound at rest to the atoms at large distance, where a frag- 
mented JSD wave function is recovered. A larger value of n 
certainly leads to a lower value of the total energy, but may im- 
prove much more the atomic energies, rather than the bonding. 

Clearly, whenever n — n* the atomization energy has to be 
referenced to the JSD calculation, even when better energies 
are provided by the RVB for the atoms. (3 It is important to 
remark here that, upon stretching the molecule to the atomiza- 
tion limit, some symmetries of the atomic wave function are 
not recovered. For instance, for a diatomic molecule such as 
iron, the total angular momentum is conserved only for the ro- 
tations around the molecular axis, so that the total angular mo- 
mentum of the atomic fragments could not be definite within 
this ansatz. Therefore, by assuming that this error does not af- 
fect the molecular bond, a correction to the atomic reference is 
necessary, given by the energy difference between the hybrid 
SD state reached upon streaching the RVB wave function and 
the fully symmetric JSD atomic state computed with the same 
primitive basis as the one used for the molecule. 

Our RVB wave function is the input for QMC simula- 
tions. We start from its optimization, and then perform vari- 
ational Monte Carlo (VMC) and diffusion Monte Carlo (DMC) 
simulations. J20I1 the latter within its recent lattice regularized 
implementation (LRDMQ. 1I21I1 The optimization method used 
here is based on the calculations of the Hamiltonian matrix el- 
ements in the space spanned by the wave function i22[ 123m 1 . re- 
cently developed to perform a constrained energy minimization 
with fixed n.[ 18] 



3. Results 

We compute the properties of the 7 A„, 7 E~ 9 E~ states for the 
neutral dimer, with the method described in Sec. 12 namely with 
the idea that only a well controlled dissociation limit can lead 
to reliable predictions for the energetics in the bonding region. 
We also compute the 8 S t 7 and 8 A g states of the anion to make a 
direct comparison with the PES. 
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Table 1: Dissociation limits of various Fe2 states for the RVB wave 
function. S D (occupation 4s 2 3d 6 ), and 5 F (occupation 4s l 3d 7 ) refer to 
neutral atom states, while 4 F (occupation 4s 2 3d 7 ) is the ground state 
of the anion. The states denoted with <2S+1) [L Z ] are non-definite angu- 
lar momentum states for the presence of other components. They are 
eigenstate of L- but not of L~. In particular 4 [0] and 5 [0] indicate the 
single occupation of the d orbitals with /- = ±2 and L = 0, 5 [1] indi- 
cates the single occupation of the orbitals with L = ±2 and L = — 1, 
while the remaining d orbitals are doubly occupied. 



Fe 2 state 


from our wf 


exact 




5 £> + 5 [l] 


5 D + 5 F 




5 D + 5 [0] 


5 D + S F 


7 A„ -» 


5 D + 5 [0] 


5 D + 5 D 


Ll g —r 


5 D + 4 [0] 


5 D + 4 F 


-» 


5 D + 4 [0] 


5 D + 4 F 



For all our calculations, we employ a neon-core pseudopo- 
tential, in order to avoid the chemically inert core electrons of 
the iron atom, and speed up the QMC simulations. We choose 
the Dolg's pseudopotential,j24j which has previously proven to 
be reliable at least for atomic QMC calculations. [25] We use a 
[8s5p6d/2slpld] contracted Gaussian basis set, which leads to 
a space spanned by 4 cr, 4 cr*, 4 n, 4 n*, 2 6, and 2 6* MO's, 
where we need to accommodate 32 electrons for the neutral 
dimer. Our primitive basis set is quite compact. However, we 
double checked the convergence in the energy differences for 
atomic calculations by extending the basis set up to 8s5p6d3f. 
A much smaller basis was used for the Jastrow factor, because 
this allows for a more efficient energy optimization. On the 
other hand, the essentially exact contribution of Jastrow-type 
dynamical correlations, which do not change the phases of the 
wave function, can be very accurately obtained with DMC or 
LRDMC. The (LR)DMC approach can be seen as a stochas- 
tic optimization of a much more general Jastrow factor which 
keeps fixed the nodes of the RVB wave function. 

For the 9 lT g , 8 S M , and & A g states n* = JV T is such that the AGP 
is an HF Slater determinant, while for the 7 A„ and the 7 E~ states 
ri* = Nf+l, and the role of non-dynamical correlations becomes 
crucial in our optimization of the AGR In fact the RVB provides 
in these latter cases a remarkable energy gain of more than 1 eV 
as compared with a more simple but much less accurate JSD 
calculation, at least as far as the agreement with experiments is 
concerned. 

All the lowest possible configurations corresponding to states 
obtained after the dissociation of our wave function are reported 
in Table [TJ It is apparent that non-dynamical correlations are 
very important in the compound whenever its total spin S is 
less than the total maximum spin of the fragments. As antic- 
ipated, a controlled atomic limit can be obtained also in this 
case, with the caveat that the total angular momentum L, around 
the molecular axis is conserved in the atomization process. This 
implies that, within our RVB wave function, the total angular 
momentum L of the fragments may not be definite in the atom- 
ization limit, and the corresponding JSD atomic reference en- 
ergies depend explicitly on the angular momentum projection 
quantum number L z . In this way, we need to perform atomic 



Table 2: The 5 F, and 5 D energies for the neutral atom, and the 4 F 
energy for the anion are reported in Hartree. The non-definite angular 
momentum JSD states are denoted with the conventions reported in 
Table [TJ They will be useful for the energy correction based on the 
dissociation limit. The variational wave function used here is a JSD 
like. From the total energies of the spin-definite states we calculate 
the 5 D — > 5 F excitations and the electron affinity, expressed in eV and 
compared with the experiment. 





LRDMC 


exp. 


S D 


-123.7819(11) 




5 F 


-123.7520(11) 




5 [0] 


-123.73986(72) 




5 [1] 


-123.71856(75) 




4 F 


-123.77731(94) 




4 [0] 


-123.76544(81) 




5 D -> 5 F (eV) 


0.81(4) 


0.87" 


4 F -> 5 D (eV) 


-0.12(4) 


0.15 a 



a From Ref. 01 



calculations for various L z and choose the appropriate ones for 
the reference to the total energy of the dimer. The results of 
our atomic calculations are reported in Table [2] The S D — > 5 F 
transition is found to be very close to the experimental value, 
while the electron affinity is off by 0.3 eV. This is largely due to 
a lack of correlation energy in the LRDMC calculations for the 
anion, since MRCI calculations done with MOLPRO[26] are in 
agreement with the experimental value. 

After performing LRDMC simulations for the iron dimer at 
different interatomic distances (going from 3.5 to 8) and sym- 
metry states, we found the results reported in Table [3] The 
vibrational frequencies of 9 I. g and & Y, g are in good agreement 
with the Leopold's experimental data. Indeed, our best QMC 
estimate for the vibrational frequency of the Fe2 ground state is 
o) e = 301(15)cm _I , which matches perfectly the value 300(15) 
cm -1 coming from PES, HI and the value 299.6 cm' 1 pro- 
vided by Raman specfroscopv. ll28ll Also for the anion dimer 
the LRDMC vibrational frequency (a> e = 210(20)cto~') agrees 
with the experimental value of 250(25) cm -1 yielded by PES. til 
Notice that our calculations correctly reproduce the softening of 
the vibrational mode going from 9 Z ? to 8 2 ? . Therefore, from 
these results we can confirm the symmetry of the peaks seen 
in the PES. Indeed, the 7 A„ state has a much higher vibrational 
frequency, incompatible with the experiment. It is interesting to 
highlight that the vibrational frequency for 7 A„ computed with 
QMC simulations agrees with those calculated by DFT meth- 
ods for the same states. Jfjl 0] To check whether the 9 £" is the 
true ground state, let us apply the correction to the energy levels 
based on the atomization limits. To do so, we add to the total 
energy of the dimer the atomic energy differences between the 
right (the exact fully symmetric JSD limit) and the middle col- 
umn (the asymptotic JSD limit that is possible to reach within 
our wave function) of Table [TJ 

The level ordering that we find is reported in Table |4] It 
turns out that after the correction the 7 A„ is the actual ground 
state, while the energy split between the 9 I. g and the 7 S^ states 
is 0.64(7) eV, in quite good agreement with the experimental 
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Table 3: The LRDMC results are reported for some states of the neu- 
tral iron dimer, and for the anion states 8 E~ and 8 A S . We calculated 
total energies at the minimum of the interatomic potential, equilibrium 
distances R e , and vibrational frequencies a> e . Some available experi- 
mental values are also reported. 





Energy 


Re 




exp aj e 




(Hartree) 


(a.u.) 


(cmr 1 ) 


(cm' 1 ) 




-247.5036(20) 


4.081(18) 


327(15) 


300(15) b 


% 


-247.5486(20) 


4.093(19) 


301(15) 


299.6 a 


7 A„ 


-247.5351(30) 


3.894(18) 


373(32) 




8 A 


-247.5585(30) 


3.908(14) 


354(24) 




u 


-247.5706(42) 


4.276(28) 


210(20) 


250(20) b 



a From Ref. IH 



b From Ref. d 



findings (0.534(4) eV).|H The correction does not change the 
ordering for the anion, and so its ground state remains the 8 £~. 
This is a quite interesting effect of the correlation, which acts on 
the d orbitals in a way which depends on the global symmetry 
and total charge of the system. This is apparent already at the 
atomic level, where the occupation of the d orbitals changes by 
going from the 4 F anion to the neutral 5 D ground states. The 
energy difference between the anion 8 £~ state and the neutral 
9 2; is 0.59(12) eV, which should be further shifted by 0.3 eV 
due to the electron affinity correction in the atomic calculations, 
because for this quantity the JSD is rather poor, as we have pre- 
viously shown in Table[2] This will lead to a difference of 0.9 
eV between the two states, again in agreement with the experi- 
mental value of 0.902(8) eV-QlJ] This suggests that our approach 
can be further improved by replacing the reference JSD atom- 
ization energies with exact atomic energies, readily available 
in quantum chemistry databases. This approach is clearly use- 
ful and practical when one needs to consider total energy dif- 
ferences between electronic states with different particle num- 
ber, like the electron affinity and ionization energies. We would 
like to stress that, although the correction based on the asymp- 
totic limits is "approximate", the accuracy level reachable in 
this way (« 0.1, maximum 0.2 eV) is below the final energy 
differences J 1811 Therefore, the ordering we propose here can 
be taken with a good confidence. 

Table 4: The results for the energy minimum are reported for the same 
states as in Table[3] but with the correction described in the text which 
takes into account the atomic limit of the wave functions. In the last 
two columns, we report the energy differences with respect to the state 
9 £~ after the correction, and the experimental values taken from the 
Leopold's experiment. 





Energy 


Corrected 


Difference 


Exp 




(Hartree) 


(Hartree) 


(eV) 


(eV) 




-247.5036(20) 


-247.5370(24) 


+0.64(7) 


+0.534(4) a 


% 


-247.5486(20) 


-247.5608(24) 


0.0 




7 A„ 


-247.5351(30) 


-247.5771(33) 


-0.44(9) 






-247.5585(30) 


-247.5703(33) 


-0.26(9) 






-247.5706(42) 


-247.5824(44) 


-0.59(12) 


-0.902(8) a 



a From Ref. |1] 



Notice that the determination of the 7 A„ as the ground state 
is also supported by the equilibrium bond lengths provided 
by our calculations. Indeed, the experimental bond length is 
3.82(4). lr29Tl which is very close to our findings for the 7 A„ 
(see Table [3]). On the other hand, the bond length of the an- 
ion has been measured only indirectly, since the unique avail- 
able data are taken from the PES, which revealed a variation of 
the equilibrium distance during the excitation from the anion to 
the neutral iron dimer. A harmonic Franck-Condon analysis of 
the vibronic band intensity profile yielded a bond elongation of 
0.15(4) a.u. on electron attachment. Ill] Now, the difference be- 
tween the 9 2 y and 8 £~ bond length of our LRDMC calculation 
amounts to 0.18 a.u., which is in perfect agreement with the 
elongation of the anion dimer measured in the PES. Therefore, 
we conclude that the states seen in the Leopold's photoelectron 
spectrum are the anion 8 £~, and the neutral and 7 Z ? . The 
ground state of the neutral dimer is however the 7 A„ state, not 
seen in the photoelectron experiment since the transition from 
a 8 E~ (the anion symmetry of the prepared initial state) to a 
7 A M symmetry is a second-order process, and so of negligible 
rate with respect to the 9 T. g and 7 E~ states, connected with the 
8 £~ by a direct photodetachment of the electron living in the 
4cr*(4s) orbital. 0] 



4. Conclusions 

We have shown that by using a cheap and nevertheless very 
accurate realization of the RVB wave function based on the MO 
expansion of the AGP part, it is possible to tackle highly de- 
bated and challenging transition metal compounds. The solu- 
tion of the iron dimer puzzle appears at end, and we strongly 
believe that many other problems - where the electron correla- 
ton plays a strong role - could be finally understood within this 
framework. 
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